La cartographie avec R

Xavier Timbeau

Pourquoi des cartes avec R ?

on peut traiter les données spatiales et réaliser des cartes avec QGIS ou d’autres logiciels

• interface adaptée, visualisation rapide

• nombreux outils de traitements

• ressources nombreuses (tuto, doc, plugin)

  • QGIS : gratuit

  • Arcgis : standard, payant, puissant (relié avec des bases de données spatiales)

  • Articque Cartes et Données payant, plus simple

  • Excel et sans doute d’autres

Pourquoi des cartes avec R ?

Avec R, on va pouvoir :

• manipuler les données spatiales comme d’autres données

  • analyse statistique, économétrie spatiale

• utiliser des outils avancés

  • visualisation, visualisation interactive

  • production de sorties pour impression, pour publication interactive (HTML5)

• automatisation des process

  • de la donnée au résultat

  • production de cartes complexes (facet, inset)

Les données spatiales (I)

des données (des lignes d’observation) plus une information géographique (une géolocalisation, une forme)

• DVF : transaction, plus un point (lon/lat) ou une parcelle (un polygone localisé)

• eurostat : données par pays, région (NUTS 1 à 3), plus information géographique (contour de la région)

• OpenStreetMap : données sur le réseau de rues et route, les batiments, des lignes, des points, des polygones, des courbes, etc…

• format dit vectoriel, manipulé avec {sf} un package R qui permet d’accéder à la librairie GDAL/GEOS qui est un standard ouvert de manipulation de données spatiales vectorielles

• les données sont sous différents formats, souvent .shp, la fonction sf::st_read s’en occupe

Les données spatiales (II)

• des rasters : un grand nombre de points ou de rectangles ou de carrés de taille identique qui recouvre un territoire

  • imagerie satellite (COPERNICUS)

  • données de pollution, méteo, température issues de modèle et d’intrapolation à partir d’observations discrètes

  • des outils spécifiques pour optimiser les traitements sur ces données et leur affichage

  • format dit raster

Le système de coordonnées

La terre est plutôt sphérique, mais elle est presque plate

coordonnées quasi-sphériques : WGS84 (World Geodetic System 1984, EPSG: 4326) tient compte de l’aplatissement de la terre, normalisé pour les GPS

Les projections

pour faire des cartes qui sont sur des plans, il faut une projection

Cylindrique	Conique	Azymutale
Mercator (conforme) Peters (équivalente) Robinson	Lambert 93 (conique conforme) (EPSG: 2154) ETRS89 Lambert Azimuthal Equal-Area projection (EPSG: 3035)	Lambert équivalente azymutale Projections gnomoniques

une projection équivalente conserve les aires, recommandé pour des cartes avec plusieurs pays ; conforme les angles (les formes), recommandé pour la topographie à de petites échelles ; aphylactique et projections uniques (autres systèmes)

Des projections dites uniques

Projection de Goode pseudo cylindrique interrompue, équivalente
Cossin (1570) sinusoïdale; équivalente
Projection gnomonique

D’un système à l’autre

Dans un sf ou un raster géographique, le système de coordonnées est défini

• sinon, on le défini (data |> st_as_sf(crs=3035)) sans se tromper

st_transform(3035) passe d’un système à la grille européenne

Le choix du système permet :

• de calculer facilement les distances ou les aires (dans une projection conforme/équivalente, localement)

• de construire une grille utile pour mixer vectoriel et raster

• modifier la représentation graphique (avec des déformations plus ou moins maîtrisées)

crs=3035 est une norme européenne. La carreau 200, COPERNICUS sont diffusés dans ce système

crs=4326 est le format le plus courant depuis l’avènement des GPS (et des calculs informatisés de changement de coordonnées)

Une première carte

Code

library(sf)
library(tmap)
library(cartogram)
library(tidyverse)
library(ofce)
library(eurostat)
bbox <-st_bbox( c(xmin=-25,ymin=32,xmax=47,ymax=71.5), crs = 4326)
sf <- get_eurostat_geospatial(
  output_class = "sf",
  resolution = "60",
  nuts_level = "all") 
st_crs(sf) <- st_crs(4326)

pib <- ofce::sna_get("nama_10r_3gdp")
pop <- ofce::sna_get("nama_10r_3popgdp")
pib <- pib |> mutate(
  geo_label = label_eurostat(geo, dic="geo", fix_duplicated = TRUE) 
)

pib <- pib  |> 
  rename(gdp = values) |> 
  filter(geo!="EU27_2020", geo!="EU28") |>  
  complete(unit, geo, time, fill=list(gdp=NA))  |> 
  mutate(pays=str_sub(geo,1,2), 
         nuts3 = str_length(str_sub(geo,3,5))==2,
         nuts2 = str_length(str_sub(geo,3,5))==1,
         nuts1 = str_length(str_sub(geo,3,5))==0) %>% 
  left_join(pop  |> dplyr::select(geo, time, pop=THS), by=c("geo","time"))

pib.sf <- left_join(
  sf  |> dplyr::select(geo=NUTS_ID, LEVL_CODE),
  pib |> 
    filter(time=="2017-01-01", unit=="PPS_EU27_2020_HAB")
  , by="geo") |>
  st_transform(3035)

mb <- mapboxapi::get_static_tiles(
  location = pib.sf |> st_transform(4326),
  zoom=2, 
  style_id = "ckjka0noe1eg819qrhuu1vigs", 
  username="xtimbeau")
mb2 <- 255 - (255-mb)/3

map <- tm_shape(pib.sf |> filter(LEVL_CODE==1), bbox=bbox)+
  tm_fill(col="gdp", style="cont", palette="-RdBu", alpha=0.75)+
  tm_shape(pib.sf |> filter(LEVL_CODE==0))+
  tm_borders()

tmap_leaflet(map)

Cette carte pose de nombreux problèmes

1. la palette de couleur : suppose un point milieu, défini une zone “ok” une zone non “ok”

2. la surface colorée dépend de la surface des territoire. Or, la donnée est une donnée de richesse par tête. Donne l’impression que l’Ile de France est petite (12m) comparée à l’Irlande (5m)

Des solutions: une autre palette est possible

Code

(tm_shape(pib.sf |> filter(LEVL_CODE==1), bbox=bbox)+
  tm_fill(col="gdp", style="cont", palette="-viridis", alpha=0.75)+
  tm_shape(pib.sf |> filter(LEVL_CODE==0))+
  tm_borders() )|>
  tmap_leaflet()

Des solutions : des “dots” plutôt que des “fill”

Code

(tm_shape(pib.sf |> filter(LEVL_CODE==1), bbox=bbox)+
  tm_dots(col="gdp", size="pop", style="cont", palette="-viridis", alpha=0.75)+
  tm_shape(pib.sf |> filter(LEVL_CODE==0))+
  tm_borders() )|>
  tmap_leaflet()

On peut faire un cartogramme

Code

library(cartogram)

cartcont <- cartogram_cont(pib.sf %>% filter(LEVL_CODE==0, !is.na(pop)),
                           weight="pop")

cartdorling <- cartogram_dorling(pib.sf %>% filter(LEVL_CODE==1, !is.na(pop)), 
                                 weight="pop",
                                 m_weight=1,
                                 k=0.33)

• Continu
• Dorling

Code

(tm_shape(cartcont, bbox=bbox)+
  tm_fill(col="gdp", size="pop", style="cont", palette="-viridis", alpha=0.75)+
  tm_borders() )|>
  tmap_leaflet()

Code

(tm_shape(cartdorling, bbox=bbox)+
  tm_fill(col="gdp", style="cont", palette="-viridis", alpha=0.75)+
  tm_shape(pib.sf |> filter(LEVL_CODE==0))+
  tm_borders() )|>
  tmap_leaflet()

une échelle log sur les couleurs

les palettes de {scico}, ici batlow

Des ressources pour faire des cartes

• Des fonds de cartes

  • IGN Contours administratifs France, régions, départements, communes, IRIS, parcelles

  • Pays européens : {eurostat}

  • Monde {rworldmap}

• Des packages R pour la manipulation des données spatiales ou leur visualisation

  • {sf}/{tidyverse}, {raster}/{terra}, {stars}

  • {tmap}/{leaflet}, {ggplot}, {mapdeck}/{mapboxapi}

• Des personnes à l’OFCE

  • Paul Malliet, Maxime Parodi (*), Xavier Timbeau